#setting proxy
Sys.setenv(http_proxy="http://10.118.113.5:3128")
Sys.setenv(https_proxy="https://10.118.113.5:3128")

#setting path
setwd("/home/barupal/Documents/R/Deseq_protocol")
getwd()
options(internet.info = 0)
source("http://www.bioconductor.org/biocLite.R")
biocLite("BiocUpgrade")
biocLite( c("ShortRead","DESeq", "edgeR") )

# downloading example data


setwd("/home/barupal/Documents/R/Deseq_protocol/data")
sri = read.csv("SraRunInfo.csv", stringsAsFactors=FALSE)
keep = grep("CG8144|Untreated-",sri$LibraryName)
sri = sri[keep,]

#downlaod sample SRA files 
fs = basename(sri$download_path)
for(i in 1:nrow(sri))
  download.file(sri$download_path[i], fs[i])

#srs to fastaq converter 
stopifnot( all(file.exists(fs)) ) # assure FTP download was successful
for(f in fs) {
  cmd = paste("fastq-dump --split-3", f)
  cat(cmd,"\n")
  system(cmd) # invoke command
}
# building a refrence index 
#bowtie2-build -f Drosophila_melanogaster.BDGP5.70.dna.toplevel.fa Dme1_BDGP5_70
#asses sequence quality 

library("ShortRead")
fqQC = qa(dirPath=".", pattern=".fastq$", type="fastq") # last two hour
report(fqQC, type="html", dest="fastqQAreport")

#collapse the intial table sri
sri$LibraryName = gsub("S2_DRSC_","",sri$LibraryName) # trim label
samples = unique(sri[,c("LibraryName","LibraryLayout")])
for(i in seq_len(nrow(samples))) {
  rw = (sri$LibraryName==samples$LibraryName[i])
  if(samples$LibraryLayout[i]=="PAIRED") {
    samples$fastq1[i] = paste0(sri$Run[rw],"_1.fastq",collapse=",")
    samples$fastq2[i] = paste0(sri$Run[rw],"_2.fastq",collapse=",")
  } else {
    samples$fastq1[i] = paste0(sri$Run[rw],".fastq",collapse=",")
    samples$fastq2[i] = ""
  }
}


#adding columns 
samples$condition = "CTL"
samples$condition[grep("RNAi",samples$LibraryName)] = "KD"
samples$shortname = paste( substr(samples$condition,1,2),
                           substr(samples$LibraryLayout,1,2),
                           seq_len(nrow(samples)), sep=".")

samples
LibraryName
bowind
fastq1
fastq2


#aling the samples 
getwd()
gf = "Drosophila_melanogaster.BDGP5.70.gtf"
bowind = "Dme1_BDGP5_70"
cmd = with(samples,
           paste0("tophat -G ", gf, " -p 5 -o", " ", LibraryName," ", bowind," ",
                  fastq1," ", fastq2))
cmd
cat(cmd,"\n")
system(cmd)
#conv
#bedtools bamtobed -i reads.bam > reads.bed erting bam into bed
#
#Organize, sort and index the BAM les and create SAM files
for(i in seq_len(nrow(samples))) {
  lib = samples$LibraryName[i]
  ob = file.path(lib, "accepted_hits.bam")
  # sort by name, convert to SAM for htseq-count
  cat(paste0("samtools sort -n ",ob," ",lib,"_sn"),"\n")
  # cat(paste0("samtools view -o ",lib,"_sn.sam ",lib,"_sn.bam"),"\n")
  cat(paste0("samtools view -o ",lib," ",lib),"\n")
  # sort by position and index for IGV
  cat(paste0("samtools sort ",ob," ",lib,"_s"),"\n")
  cat(paste0("samtools index ",lib,"_s.bam"),"\n\n")
}

#count with htseq
samples$countf = paste(samples$LibraryName, "bed", sep=".")
gf = "Drosophila_melanogaster.BDGP5.70.gtf"
cmd = paste0("bedtools bamtobed -i ", samples$LibraryName,"/accepted_hits.bam",  " > ", samples$countf , "\n")
#cmd = paste0("htseq-count -s no -a 10 ", samples$LibraryName, "_sn.sam ",
#               gf," > ", samples$countf)
cmd
cat(cmd,"\n")
system(cmd)

